function TSMC_Focus_Map
%#ok<*NASGU>
%#ok<*ASGLU>
%#ok<*UNRCH>

GoldenInputImage='C:\Users\Les Frackowiak\Desktop\Tests\TSMC\Coarsely Aligned\No Letters\46-G.jpg';
SampleInputImage='C:\Users\Les Frackowiak\Desktop\Tests\TSMC\Coarsely Aligned\No Letters\46-S.jpg';
UseBackgroundEqualization=1;
BackgroundBlurSize=15;
BackgroundBlurSigma=10;
UseGrayscaleEqualization=1;
AbsoluteGrayscaleValue=150;
TrimAmount=30;
FilterType=0;
% FilterType = 0 : None
% FilterType = 1 : Median
FilterSize=11;
% Must be odd integer
InterpolationMethod=1;
% InterpolationMethod = 1 : Multiplication
FocusCapValue=2;
SmoothType=1;
% FilterType = 0 : None
% FilterType = 1 : Average
SmoothSize=50;

%%%% IMAGE LOADING %%%%
GoldenImage=imread(GoldenInputImage);
SampleImage=imread(SampleInputImage);
if size(size(GoldenImage),2)==3
    GoldenGray=rgb2gray(GoldenImage);
else
    GoldenGray=GoldenImage;
end
if size(size(SampleImage),2)==3
    SampleGray=rgb2gray(SampleImage);
else
    SampleGray=SampleImage;
end
[GoldenHeight,GoldenWidth]=size(GoldenGray);
[SampleHeight,SampleWidth]=size(SampleGray);
MinWidth=min(GoldenWidth,SampleWidth);
MinHeight=min(GoldenHeight,SampleHeight);
GoldenGray=GoldenGray(1:MinHeight,1:MinWidth);
SampleGray=SampleGray(1:MinHeight,1:MinWidth);
GoldenGray=GoldenGray(TrimAmount:(MinHeight-TrimAmount),TrimAmount:(MinWidth-TrimAmount));
SampleGray=SampleGray(TrimAmount:(MinHeight-TrimAmount),TrimAmount:(MinWidth-TrimAmount));

%%%% IMAGE PREP %%%%
if UseBackgroundEqualization
    Filter=fspecial('gaussian',BackgroundBlurSize,BackgroundBlurSigma);
    GI_Filtered=imfilter(GoldenGray,Filter,'symmetric');
    SI_Filtered=imfilter(SampleGray,Filter,'symmetric');
    GI_Filtered_Min=min(min(GI_Filtered));
    SI_Filtered_Min=min(min(SI_Filtered));
    GoldenGray=GoldenGray-(GI_Filtered-GI_Filtered_Min);
    SampleGray=SampleGray-(SI_Filtered-SI_Filtered_Min);
end
if UseGrayscaleEqualization
    GI_AverageValue=mean2(GoldenGray);
    SI_AverageValue=mean2(SampleGray);
    GoldenGray=GoldenGray*(AbsoluteGrayscaleValue/GI_AverageValue);
    SampleGray=SampleGray*(AbsoluteGrayscaleValue/SI_AverageValue);
end

%%%% 1D CONVERSION %%%%
GoldenHorizontal=std(double(GoldenGray),0,1);
SampleHorizontal=std(double(SampleGray),0,1);
GoldenVertical=std(double(GoldenGray),0,2);
SampleVertical=std(double(SampleGray),0,2);
if FilterType==1
    GoldenHorizontal=medfilt1(GoldenHorizontal,FilterSize);
    SampleHorizontal=medfilt1(SampleHorizontal,FilterSize);
    GoldenVertical=medfilt1(GoldenVertical,FilterSize);
    SampleVertical=medfilt1(SampleVertical,FilterSize);
end
if InterpolationMethod==1
    GoldenMap=GoldenVertical*GoldenHorizontal;
    SampleMap=SampleVertical*SampleHorizontal;
end
FocusMap=SampleMap./GoldenMap;
FocusMap=abs(FocusMap-1);
for i=1:MinWidth-2*TrimAmount+1
    for j=1:MinHeight-2*TrimAmount+1
        Value=FocusMap(j,i);
        if (Value>FocusCapValue)
            FocusMap(j,i)=FocusCapValue;
        end
    end
end
if SmoothType==1
    FocusMap=imfilter(FocusMap,fspecial('average',SmoothSize),'replicate');
end
figure, imagesc(FocusMap,[0 FocusCapValue])
end


